Relaxation and Glassy Dynamics in Disordered Type-II Superconductors 
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We study the non-equihbrium relaxation kinetics of interacting magnetic flux lines in disordered 
type-II superconductors at low temperatures and low magnetic fields by means of a three-dimensional 
elastic line model and Monte Carlo simulations. Investigating the vortex density and height auto- 
correlation functions as well as the flux line mean-square displacement, we observe the emergence 
of glassy dynamics, caused by the competing effects of vortex pinning due to point defects and 
long-range repulsive interactions between the flux lines. Our systematic numerical study allows us 
to carefully disentangle the associated different relaxation mechanisms, and to assess their relative 
impact on the kinetics of dilute vortex matter at low temperatures. 
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I. INTRODUCTION 

In this paper, we report an investigation of the non- 
equilibrium relaxation kinetics in the vortex glass phase 
of layered disordered type-II superconductors. Since 
Struik's original investigationsji many glassy systems 
have been found to exhibit physical aging phenomena, 
which have attracted considerable interest during the 
past decades.— Recently, it has been realized that glass- 
like relaxation and aging can in fact be found in many 
other systems.'- - Glassy materials feature extremely long 
relaxation times which facilitates the investigation of ag- 
ing phenomena in real as well as in numerical experi- 
ments. Our definition of physical aging here entails two 
fundamental properties: First, we require relaxation to- 
wards equilibrium to be very slow, typically character- 
ized by a power law decay, observable in a large acces- 
sible time window tmic 4; i <C ieq; here <mic denotes an 
appropriate short microscopic time scale, whereas teq is 
the much larger equilibration time for the macroscopic 
system under consideration. Second, a non-equilibrium 
initial state is prepared such that the kinetics is rendered 
non-stationary; thus, time-translation invariance is bro- 
ken, and two-time response and correlation functions de- 
pend on both times s and t > s independently, not just 
on the elapsed time difference t — s. In this context, s 
is often referred to as waiting time, and t as observation 
time. In addition, in the limit t^s many aging systems 
are characterized by the emergence of dynamical scaling 
behavior.— 

The physics of interacting vortex lines in disordered 
type-II superconductors is remarkably complex and has 
been a major research focus in condensed matter physics 
in the past two decades. It has been established that the 
temperature vs. magnetic-field phase diagram displays a 
variety of distinct phases.^ A thorough understanding of 
the equilibrium and transport properties of vortex matter 
is clearly required to render these materials amenable to 
optimization with respect to dissipative losses, especially 
in (desirable) high-field applications. Investigations of 
vortex phases and dynamics have in turn enriched con- 
densed matter theory, specifically the mathematical mod- 



eling and description of quantum fluids, glassy states, 
topological defects, continuous phase transitions, and dy- 
namic critical phenomena. An appealing feature of dis- 
ordered magnetic flux line systems is their straightfor- 
ward experimental realization which allows direct com- 
parison of theoretical predictions with actual measure- 
ments. The existence of glassy phases in vortex mat- 
ter is well-established theoretically and experimentally.-Sii 
The low-temperature Abrikosov lattice in pure flux line 
systems is already destroyed by weak point-like disor- 
der (such as oxygen vacancies in the cuprates). The 
first-order vortex lattice melting transition of the pure 
system^ is then replaced by a continuous transition into 
a disorder-dominated vortex glass phasefS,"^ Here, the 
vortices are collectively pinned, displaying neither trans- 
lational nor orientational long-range order. In addition, 
there is now mounting evidence for a topologically or- 
dered dislocation-free Bragg glass phase at low magnetic 
fields or for weak disorder j^'^^i^'^" — and an intriguing in- 
termediate multidomain glass state has been proposed. — 

Unambiguous signatures of aging in disordered vortex 
matter have also been identified experimentally: For ex- 
ample, Du et al. recently demonstrated that the volt- 
age response of a 2H-NbSe2 sample to a current pulse 
depended on the pulse duration^^ (see also Ref. [18|). 
Out-of-equilibrium features of vortex glass systems re- 
laxing towards their equilibrium state were studied some 
time ago by Nicodemi and Jensen through Monte Carlo 
simulations of a two-dimensional coarse-grained model 
system'^ however, this model applies to very thin 
films only since it naturally disregards the prominent 
three-dimensional flux line fluctuations. More recently, 
three-dimensional Langevin dynamics simulations of vor- 
tex matter were employed by Olson et al^ and by 
Bustingorry, Cugliandolo, and Dommguea^iiS^ (see also 
Ref. '23','24]) in order to investigate non-equilibrium re- 
laxation kinetics, with quite intriguing results and indica- 
tions of aging behavior in quantities such as the two-time 
density-density correlation function, the linear suscep- 
tibility, and the mean-square displacement. Roma and 
Domfnguez extended these studies to Monte Carlo simu- 
lations of the three-dimensional gauge glass model at the 
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critical temperature 

We remark that it is generally crucial for the analysis of 
out-of-equilibrium systems to carefully investigate alter- 
native microscopic realizations of their dynamics in order 
to probe their actual physical properties rather than ar- 
tifacts inherent in any mathematical modeling. Indeed, 
different mathematical and numerical representations of 
non-equilibrium systems rely on various underlying a pri- 
ori assumptions that can only be validated a posteriori. 
It is therefore imperative to test a variety of different 
numerical methods and compare the ensuing results in 
order to identify those properties that are generic to the 
physical system under investigation. In this paper we 
employ Metropolis Monte Carlo simulations for a three- 
dimensional interacting elastic line model to investigate 
the relaxation behavior in the physical aging regime for 
systems with uncorrelated attractive point defects. 

We strive to employ parameter values that describe 
high-Tc superconducting materials such as YBCO, and 
limit our investigations to low magnetic fields and tem- 
peratures (typically 10 K) in order for our disordered 
elastic line model to adequately represent a type-II super- 
conductor with realistic material characteristics. Thus 
we address a parameter and time regime wherein the 
slow dynamics is dominated by the gradual build-up of 
correlations induced by an intricate interplay of repulsive 
vortex interactions and attractive point pinning sites. 

Our work differs in crucial aspects from other recent 
studies (mil]. As in Refs. [iSiO], we consider the 
in type-II superconductors physically relevant situation 
where all defects serve as genuinely attractive and local- 
ized pinning sites for vortices, in the sense that they lo- 
cally reduce the chemical potential (or equivalently, sup- 
press the superconducting transition temperature). Our 
pinning potential landscape is therefore characterized by 
large flat regions in space, where the vortices feel no pin- 
ning force, interspersed with small attractive potential 
minima of extension bo, much smaller than the London 
penetration depth Aab that sets the vortex-vortex inter- 
action range. This is to be contrasted with the model 
used in Refs. [2lll2^ . which is rather motivated by studies 
of interfaces in random environments that are described 
by Gaussian distributions for the disorder strengthj^i^ 
Consequently these models inevitably incorporate both 
attractive and repulsive disorder, which can be viewed as 
mimicking a sample with a very high density of point de- 
fects. Alternatively, such a coarse-grained representation 
of pinning centers forming a continuous disorder land- 
scape certainly becomes appropriate at elevated temper- 
atures near Tc, since then the pinning range is set by the 
coherence length ^ab, which diverges as the critical point 
is approached. Thus, a random medium description is 
best-suited for investigations of critical phenomena, and 
also more easily amenable to field-theory representations. 
At low temperatures, however, where fab < bo, our mod- 
eling of the localized pinning centers appears more re- 
alistic, and we furthermore remark that in this scenario 
repulsive defects would introduce different physics in the 



non-equilibrium relaxation and aging kinetics of vortices 
in superconductors, such as flux bunching in regions de- 
void of such disorder. We therefore carefully exclude any 
repulsive pinning sites. In addition, the temperatures 
used in Refs. .21.22.] appear to be considerably higher 
than those studied in our present work. Further differ- 
ences can be found in the length of the vortex lines (in 
Refs. [21,22,] rather short lines were considered), in the 
boundary conditions, as well as in the initial preparation 
of the system. 

We characterize the aging properties of the interacting 
and pinned flux lines through several different two-time 
quantities, namely the vortex density-density autocor- 
relation function, the flux line height-height autocorre- 
lation function, and the transverse vortex mean-square 
displacement. Investigating the influence of weak point 
defects, we find that the non-equilibrium relaxation prop- 
erties of magnetic flux lines in disordered type-II super- 
conductors are governed by various crossover effects that 
reflect the competition between pinning and repulsive in- 
teractions. In the long-time limit and for not too large 
pinning strengths, the dynamics is manifestly similar to 
that observed in structural glasses. 

The structure of this paper is as follows: In Section HIl 
we describe our model and the Monte Carlo simulation 
algorithm, and deflne the quantities of interest for our 
study. Our data and principal results are presented in 
Section IIIII In order to disentangle the different con- 
tributions to the non-equilibrium relaxation dynamics of 
our system, we first separately elucidate the effects of at- 
tractive pinning centers and of long-range vortex-vortex 
interactions, before we endeavour to analyze and under- 
stand their intriguing interplay as reflected in the vortex 
system's relaxation kinetics. Finally, we discuss our flnd- 
ings in Section HVl and compare them with other studies. 



II. MODEL AND SIMULATION PROCEDURE 



A. Effective model Hamiltonian 



We consider a three-dimensional vortex system in the 
London limit, where the London penetration depth is 
much larger than the coherence length. We model the 
vortex motion by means of an elastic flux line free en- 
ergy described in Ref. [131, see also, e.g., Refs. j28i - t32| . 
The system is composed of N flux lines in a sample of 
thickness L. The effective model Hamiltonian is de- 
fined in terms of the flux line trajectories Yj{z), with 
j = 1, . . . ,N, and consists of three components, namely 
the elastic line tension, the repulsive vortex- vortex inter- 
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action, and the disorder-induced pinning potential: 




Here, the elastic line stiffness is ei « F^^eo ln(Aab/Cab), 
with Aab and fab denoting the London penetration depth 
and coherence length in the crystallographic ab plane 
(we assume the magnetic field along the c axis), and 
the anisotropy parameter (effective mass ratio) F^^ ~ 
M^^/Mc- The energy scale is set by eo = (0o/47rAab)^, 
where 00 = hc/2e is the magnetic flux quantum. The 
expression for the elastic line energy in Eq. ([T]) is valid 
in the limit \dYj{z)/dz\'^ ^ F^. The purely in-plane 
repulsive interaction potential (consistent with the ex- 
treme London limit) between flux line elements is given 
by the modified Bessel function of zeroth order, V{r) = 
2eo-ft'o(''/Aab), which diverges logarithmically as r — )■ 0, 
and decreases exponentially for r ^ Aab- These vortex 
interactions are truncated at half the system size, which 
is in turn chosen sufficiently big such that numerical arti- 
facts due to this cut-off length are minimized. We model 
point pinning centers through square potential wells with 
radius &o and strength C/q at No defect positions. (For 
additional details, see Refs. Plsaj.) 

B. Numerical parameter values 

Our simulation parameter values^ were chosen cor- 
responding to typical material parameters for YBCO as 
listed in appendix D of Ref. [l^] . In the following, lengths 
and energies are reported relative to the effective defect 
radius feo and interaction energy scale eo (using cgs units) , 
and time in Monte Carlo steps (MCS), where one MCS 
correspond to NL proposed updates of the flux line ele- 
ments, with N the number of flux lines and L the num- 
ber of layers. We set the pinning center radius 6o = 35 
A, anisotropy parameter F^^ = 1/5, and, as is appro- 
priate at low temperatures, Aab = 346o ~ 1200 A, and 
^ab = 0.36o « 10.5 A. Then eo « 1.9 x 10"^ (in cgs units 
of energy / length) , and the energy scale in the line ten- 
sion term becomes ei w O.lSeo. We systematically vary 
the pinning strength Uq between and 0.2eo. Usually, 
our simulations are performed at temperature T = 10 K, 
which corresponds to kBT/eobo « 0.002. Thermal exci- 
tation energies are thus small compared to the elastic and 
pinning energies, and at equilibrium we therefore expect 
the system to be deep in the glassy regime. We do not 
allow for flux line cutting and reconnection processes in 
our simulations of a low-temperature and dilute vortex 
system. 



randomly placed measure C(t,s) 

vortices 

^ T=10K T=10K ^ 
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FIG. 1: Sketch of the measurement protocol. At t = 0, 
straight vortex lines are initialized far out of equilibrium by 
placing them randomly in the system at T = 10 K in the pres- 
ence of weak point defects. The vortex lines are then allowed 
to relax for different waiting times s, before various two-time 
quantities, such as the height autocorrelation function C{t, s), 
are measured. 



C. System preparation and simulation protocol 



We apply the standard Metropolis Monte Carlo sim- 
ulation algorithm in three dimensions with a discretized 
version of the above effective Hamiltonian ([T])f2^ The sys- 
tem contains = 16 flux lines in L layers, with a dis- 
tance bo between consecutive layers; and an equal number 
Njj/L = 1116 of point pinning centers that are however 
randomly distributed within the layer, with mean sepa- 
ration ~ 96o; in comparison, the triangular vortex lat- 
tice spacing would be 78.56o in our dilute system. We 
apply periodic boundary conditions in all three space di- 
rections, as we are mainly interested in bulk properties. 
This is to be contrasted with Refs. [2l1l22j . where free 
boundary conditions were used along the c axis. We have 
systematically changed L between 10 and 2560 in order 
to carefully monitor finite-size effects. The in-plane sys- 
tem size is [Lx,Ly] = [-^ x 8Aab,8Aab]; the dimensions 
of the xy plane were chosen such that in the absence of 
disorder the system accommodates a regular triangular 
flux lattice. In the absence of defects, we have tested that 
initially randomly placed vortices properly equilibrate to 
form a triangular Abrikosov flux lattice. We have also 
checked that there are no appreciable effects due to the 
sharp cut-off of the vortex interactions at 4 Aab = Ly/2. 

In order to investigate aging phenomena in the system 
with uncorrelated point disorder, the vortices are pre- 
pared in an out-of-equilibrium state: Straight flux lines 
are initially (at t = 0) placed at random locations in 
the system. The vortex lines are subsequently allowed 
to relax at the temperature T = 10 K for a duration s, 
the 'waiting' time, before we start measuring two-time 
quantities for t > s, see Fig. [T] (This is again different 
to Refs. [2lH22j . where the vortex lines were equilibrated 
at high temperatures inside the vortex liquid phase be- 
fore the subsequent quench to lower temperatures.) Our 
waiting times extend up to s = 51200 MCS, whereas the 
total length of a simulation run is typically t = 10s. 
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D. Measured quantities 

Aging phenomena can generally be adequately char- 
acterized through the study of two-time quantities. In 
our work we put special emphasis on a range of observ- 
ables that allow us to rather comprehensively monitor 
the distinct relaxation processes in vortex matter that 
originate from pinning to attractive point defects and re- 
pulsive interaction forces, respectively, and their intricate 
competitive interplay. 

The height-height autocorrelation function and mean- 
square displacement represent two quantities that are 
routinely studied in the context of interface fluctuations 
and non-equilibrium growth processesi^^— i^"— Sepa- 
rating the time-dependent position of the flux line j in 
the zth layer into its x and y components, rj{z,t) = 
(^Xj{z,t),yj{z,t)) , the two-time height- height autocorre- 
lation function can be written as 

L N 



z=lj=l 



L N 



LN 

L N 

where x{t) — ^ ^ Xj{z,t), and similarly for y{t). 

The brackets (• • • ) here denote both an average over the 
noise history, i.e., over the sequential realizations of ran- 
dom number sequences, as well as a configurational aver- 
age over defect distributions and initial positions of the 
straight vortex lines at the outset of the simulation runs. 
The two-time mean-square displacement in the xy planes, 
transverse to the external magnetic field, can similarly be 
cast in the form 

L N 
z=l 3 = 1 

+ (^{y,{z,t)-y,{z,s)y)] . (3) 



We remark that other related quantities that essentially 
contain the same information, are the two-time roughness 
function or the two-time structure factor.— 

Unfortunately, both the height autocorrelation and the 
mean-square displacement are probably not easily acces- 
sible in experiments on type-II superconductors, except 
perhaps through low-angle neutron scattering. Much 
better suited for an experimental study is likely the (con- 
nected) two-time vortex density-density autocorrelation 
function that can formally be written as 



(4) 



where p{r, t) represents the local flux density per unit 
area at position r, with constant uniform average 
(p(r,<)^ = p. Following an initially random placement. 
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FIG. 2: Sketch (two-dimensional cross section) of vortex line 
elements (small solid dots) and the associated circles with 
radius abo (large open circles) at different observation times.— 
(a) At t = s, each vortex line element by construction resides 
inside its own circle and gives a count rii (t = s) = Ns 
(= 10 here), (b) At later times t > s, the repulsive vortex 
interactions cause the flux line elements to move away from 
their initial positions. This results in a smaller occupation 
number '^^ni{t > s) < Ns {~ 5 in this example), (c) At 
long times f 2> s, it is possible that all vortex line elements 
have left the circle, ni{t ^ s) — 0, which results in a complete 
decorrelation. For instance, the quantity ni{s)ni{t) in the 
pictures at these three times is evaluated to be 10 at f = s, 5 
at t > s, and for t ^ s. 



the repulsive vortex interactions cause positional rear- 
rangements, such that one would expect a temporal de- 
cay of the density autocorrelation function. In our sim- 
ulations, we realize the vortex density autocorrelation 
function in the following wayj^ As before we start with 
randomly placed straight vortex lines at t = and let the 
system subsequently relax up to waiting time s. A den- 
sity count for the vortex line elements is then generated 
by setting a circular area, with a radius equal to a bo at 
the location of each vortex line element i = (j, z)att = s. 
Typical values for a range from 0.05 to 0.20. As time t 
elapses, we count the number of vortex line elements still 
in their circles, generating a time sequence of occupation 
numbers ni{t), with = or 1, and ni(t = s) = 1 by 
construction. Due to the repulsive vortex interactions, 
flux line elements tend to move away from their initial 
positions, whence ni(t > s) can be at a later time if 
the vortex line clement leaves the prescribed circle. In 
the presence of pinning centers, vortex line elements will 
become trapped inside the defects over a long time, caus- 
ing ni(t > s) to preferentially remain 1. This quantity 
is then averaged over the Ns = LN different vortex seg- 
ments and many distinct defect distributions and initial 
configurations, yielding 



1 



(5) 



Fig. [2] illustrates the algorithm for calculating the density 
autocorrelation function.— We have checked the results 
for different values of a and found that within a rea- 
sonable range the precise choice of a does not affect the 
results in the long-time aging regime where t ^ s. 



5 



III. RELAXATION PROCESSES 

In order to fully understand the non-equilibrium re- 
laxation processes and aging phenomena in disordered 
type-II superconductors at low temperatures, we found 
it imperative to carefully disentagle the dynamical con- 
tributions originating from the repulsive interactions be- 
tween the vortex lines and from their pinning to attrac- 
tive point defects. We start our discussion with free flux 
lines, mainly in order to validate our code by comparing 
our data with the theoretically expected behavior and 
earlier work. We will then separately consider the effects 
of attractive point pins and of the long-range repulsive 
vortex interactions, before we at last venture to study the 
interplay of these two competing mechanisms to induce 
or relax correlations in the system. 



A. Free elastic line 

The relaxation kinetics of a single free elastic vortex 
line constitutes a valuable benchmark to check our pro- 
gram as this case can be easily understood by recalling 
that in the presence of thermal noise a fluctuating in- 
terface that tries to minimize its line tension should be 
described in the continuum limit by the linear Edwards- 
Wilkinson equation. '^^ As the fluctuations in the trans- 
verse X and y directions are independent random vari- 
ables for our free line, we expect the results for the free 
vortex to be described by the one-dimensional version of 
that well-known stochastic equation {h below stands for 
either x or y): 



dh{z,t) 
dt 



vdlh{z,t)+T]{z,t) , 



(6) 



where 77(2;, t) represents a Gaussian white noise with zero 
mean and covariance {ri{z,t)ri{z' ,t')) = ^S{t — t')6{z — 
z'), V is the line stiffness (equal to 1\ here), and T the 
temperature of the heat bath. The Edwards- Wilkinson 
equation, as well as a range of microscopic models belong- 
ing to the same dynamic universality class, has been stud- 
ied extensively. Starting from a straight line, one first ob- 
serves a short-time regime with uncorrelated fluctuations, 
which is rapidly replaced by a correlated intermediate- 
time interval characterized by a non-trivial power law 
increase of the line roughness. After a crossover time 
that algebraically depends on the system size, this corre- 
lated regime flnally reaches the steady-state or saturation 
regime. 

Two-time quantities have also been studied in the con- 
text of the Edwards- Wilkinson equation^^i^ii^, facili- 
tated by the fact that a full analytical analysis is pos- 
sible for the linear stochastic equation ([6]). For ex- 
ample, in the correlated regime of the one-dimensional 
Edwards- Wilkinson equation the following exact expres- 
sion for the height-height autocorrelation function has 
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FIG. 3: (Color online) Various two-time quantities for the 
free elastic vortex line obtained in systems with (a-c) L = 10 
and (d-f) L = 2560 planes, and averaged over typically 200 
independent Monte Carlo simulation runs; (a) and (d): mean- 
square displacement, (b) and (e): height-height autocorrela- 
tion function, (c) and (f): vortex density-density autocorre- 
lation function (measured with a — 0.05). Data obtained 
for different waiting times are shown. For L = 10 layers the 
system rapidly reaches the steady state, and time-translation 
invariance is recovered, i.e., the two-time quantities only de- 
pend on the time difference t — s. For L = 2560 aging and 
dynamic scaling prevail throughout the simulation time win- 
dow. The full line in (e) indicates the exact expression ((Tj 
derived from the Edwards-Wilkinson equation. The dashed 
lines in (d) and (f) indicate the predicted asymptotic power 
laws with exponents 1/2 and —1/2, respectively. Here and 
in the following figures error bars are much smaller than the 
symbol sizes. 



been derived)^ 



1/2 



1/2 



(7) 



where Cq is a known constant. The detailed crossover 
properties of two-time quantities in the region between 
the correlated and saturated regimes have been carefully 
investigated in Ref. j23| . 

In Fig. [3l we display our Monte Carlo simulation re- 
sults for our elastic vortex line model when both the vor- 
tex interaction and the defect pinning are switched off, 
i.e., only the flrst contribution in ^ is retained. One 



immediately notices a striking difference between the be- 
havior of a "thin film" composed of only a few layers 
(such as L = 10, see Figs.|3K-|3h) and "bulk" systems con- 
sisting of many layers (L — 2560, see Figs.[3Ji-15f). In the 
former case, the system rapidly evolves into the steady- 
state regime, yielding two-time quantities that only de- 
pend on the elapsed time difference t — s. As a result, 
the transverse displacements in the x and y directions 
perform simple random walks, as is revealed by the lin- 
ear increase of the mean-square displacement with time, 
see Fig. [5^. For the larger bulk system, the correlated 
regime persists throughout the duration of our simula- 
tions, and both waiting and observation times reside well 
within that extended intermediate regime. This gives rise 
to aging and dynamical scaling: Time-translation invari- 
ance is broken, and all the two-time quantities display 
full-aging scaling.^ For each two-time observable we find 
the following scaling behavior (here given for the height 
autocorrelation function C(t,s)): 

C{t,s)^s-'fc{t/s) , (8) 

where b represents an aging scaling exponent and fciu) 
denotes an associated scaling function that follows a 
power law decay for large arguments. For the height- 
height autocorrelation we have b — — l/2i^ In Fig. [SJ? 
we explicitly compare our numerically determined scal- 
ing function with the expression ([7]) resulting from the 
direct solution of the Edwards- Wilkinson equation (full 
line), and obtain perfect agreement. 

Summarizing, we see that the free vortex line fiuctua- 
tions are indeed aptly described by the one-dimensional 
Edwards- Wilkinson equation (jG)). We also observe a 
strong dependence on the system's extension L in the 
magnetic field direction, i.e., the vortex length: On the 
time scale of our simulations, the stationary regime is al- 
most immediately reached when the system consists of 
only a few layers; in contrast, for larger bulk systems 
aging and dynamical scaling are observed easily. This 
points to quite distinct relaxation behavior in thin su- 
perconducting films and thicker bulk samples. We de- 
cided to avoid the additional complications stemming 
from the crossover between the correlated and steady- 
state regimes in our present study, and to rather focus 
on system sizes sufficiently large that no finite-size ef- 
fects (no crossover to the steady state) are observed on 
the accessed time scales. Properties of smaller systems 
and possible experimental consequences for thin super- 
conducting films will be discussed in a separate publica- 
tion. 



B. Pinning without interactions 

Intuitively, one anticipates pinning centers to strongly 
influence the thermal fluctuations of our elastic flux lines. 
Indeed, attractive forces emanating from the pinning cen- 
ters will tend to localize vortex segments, and thus ulti- 
mately suppress thermal fluctuations. Depending on the 
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FIG. 4: (Color online) Mean-square displacement B{t, 0) vs. 
time t (in MCS) for different values of the pinning strength 
p, for systems of size L — 640. In the initial time regime, the 
pinning centers attract the vortex segments as revealed by an 
increase of the slope of B{t, 0) as compared to a pure system. 
At later times the localization of the flux line elements induced 
by the pinning yields a strong decrease of this slope. The data 
shown result from averaging over typically 100 independent 
runs. The dashed lines indicate the times t = 100, 400, 1600, 
and 6400 MCS, see Fig. below. 

pinning strength, flux line elements will end up spending 
an appreciable amount of time close to a pinning cen- 
ter. Therefore, compared to freely fluctuating lines, a 
marked increase of correlations as function of time must 
be expected. 

Before we proceed to analyze the influence of pinning 
centers in more detail, we need to stress that we ex- 
clusively consider attractive point defects, in accordance 
with the physics of disordered type-II superconductors 
in the low-temperature regime. A recent study [2^ ad- 
dressed the relaxation and aging properties of elastic 
lines subjected to a random potential, corresponding to 
both attractive and repulsive pinning centers. Whereas 
a Gaussian disorder strength distribution certainly is a 
good model for disordered ferromagnets, its relevance for 
relaxation processes in disordered type-II superconduc- 
tors at low temperatures is less obvious. 

Let us start by looking at the mean-square displace- 
ment B{t,0), with s — 0, which gives a measure of 
the (squared) distance traveled by the flux line elements 
since the initial preparation of the system. In Fig. |4] we 
compare the behavior of a free elastic line with that of 
flux lines subject to pinning centers of various strengths 
p = Uq/^o- The presence of attractive pins clearly gives 
rise to different regimes. The flux lines are rapidly at- 
tracted by the point defects, which yields an increase 
of the slope in the log-log plots of B{t, 0) vs. time t. 
This continues until some pinning strength-dependent 
crossover time at which the slope decreases even below 
the value of the free line, signifying the confinement of lo- 
calized vortex segments to the vicinity of the pins. As one 
would expect, this crossover time decreases for increasing 
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pinning strengths. For p > 0.05, B{t,0) remains essen- 
tially unchanged, which indicates that for non-interacting 
lines there exists a critical pinning strength above which 
thermal fluctuations are no longer sufficiently strong to 
allow the vortex line elements to escape from the defects. 

These different regimes also manifest themselves when 
two-time quantities are considered, as seen in Fig. [SJ 
where we have plotted the data according to the free-line 
scaling behavior. Of course, it is not to be expected that 
these scaling laws remain valid when attractive defects 
are added to the system, but this representation of our 
data facilitates the following discussion. We first remark, 
see Figs. [S^ and [SJl, that the change of the slope oiB{t, 0) 
translates into deviations of the mean-square displace- 
ment B{t, s) from the free-line scaling that can be readily 
understood. For example, for p = 0.01 the time intervals 
[25,250], [100,1000], and [400,4000], used to compute 
B{t, s) for the waiting times s — 25, 100, and 400, re- 
spectively, correspond to the time regime with increasing 
local slopes of B{t, 0), compare Fig.lH This yields a shift 
of ln{s^'^'^ B{t, s)) to higher values. As the crossover time 
of B{t,0) lies inside the interval [1600,16000], the con- 
verse behavior is observed for s = 1600 and even larger 
waiting times, with a shift of \n{s~'^^^B{t, s)) to lower 
values. This effect is more pronounced for larger pinning 
strengths, since then the crossover of B{t, 0) takes place 
earlier. 

The strongest influence of point defects and largest de- 
viations from free elastic lines are observed in the height 
autocorrelation function, Figs.[5}D and[5^. As the flux line 
elements are trapped by the pinning centers, their trans- 
verse in-plane displacements become diminished, which 
leads to an increase of the correlations as a function of 
waiting time. In addition, the decay of C{t, s) as a func- 
tion of t is much slower for larger values of s. For larger 
pinning strengths and long waiting times we even observe 
non-monotonic temporal evolution, as the trapped fiux 
lines experience an increase of the height correlations. 

Finally, the vortex density-density autocorrelation 
turns out to be the least sensitive among our two-time 
observables to the presence of pinning centers, at least 
for comparatively small waiting times, see Figs. [SJ: and 
[Sf. Indeed, for moderate values of s one still observes the 
free-line scaling behavior; this is of course a consequence 
of our prescription for the computation of this correlation 
function, namely setting a circular area with fixed radius 
around every vortex line element at time s: As long as 
only a few line elements are captured by a point defect, 
the scaling of Cy{t, s) remains approximately unchanged. 
Only when the majority of the vortex segments become 
trapped, does this localization induce a strong enhance- 
ment of the vortex density-density correlations. 

To conclude this section, we note that the non- 
equilibrium relaxation physics is drastically different 
when both attractive and repulsive pinning centers are 
implemented. As studied in Ref. [l^j (see also Refs. [H, 
Isoj ) an elastic line in a random potential is characterized 
by a time-dependent correlation length that crosses over 
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FIG. 5: (Color online) Various two-time quantities for non- 
interacting vortex lines subject to attractive point defects of 
strengths (a-c) p = 0.01 and (d-f) p = 0.05: (a) and (d) 
mean-square displacement; (b) and (e) height-height auto- 
correlation function; (c) and (f) vortex density-density au- 
tocorrelation (with Q = 0.05). The shown data result from 
averaging over typically 100 independent runs in systems of 
size L = 640. For direct comparison with Fig. [3l the data 
are plotted according to the scaling properties of free elas- 
tic lines. Whereas an approximate scaling prevails for small 
waiting times (especially for the mean-square displacement 
and the vortex density-density autocorrelation), strong devi- 
ations emerge for larger waiting times. In (e) the apparent 
collapse of the height autocorrelation function for the largest 
waiting times is merely caused by the scale used in that figure. 
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FIG. 6: (Color online) (a) Mean-square displacement and 
(b) height-height autocorrelation function for non-interacting 
vortex lines subject to both attractive and repulsive point 
defects with strengths drawn from the square distribution 
[—0.01, 0.01]. The shown data result from averaging over typ- 
ically 100 independent runs in systems of size L = 1280. The 
data display simple aging scaling, with different exponents for 
the two quantities. 



from an early time power law growth to an asymptotic 
logarithmic growth. Consequently, two-time quantities 
display an apparent simple aging scaling with effective 
exponents that depend on temperature and on the ran- 
domness. We have verified that we obtain similar results 
as in Ref. when using both attractive and repulsive 
defects in our model and Monte Carlo algorithm. In- 
deed, as shown in Fig. [Bl the time-dependent correlation 
length gives rise to simple aging scaling of our two-time 
quantities, with effective exponents that display a de- 
pendence on temperature and on the distribution of the 
pinning strengths. These crossover features also capture 
most of the relevant properties of disordered ferromag- 
nets undergoing phase ordering42ri^ In Refs. 2lll2^ . dis- 
ordered type-II superconductors in the low-temperature 
phase were modeled by a corresponding model with ran- 
dom pins that are either attractive or repulsive. How- 
ever, the physical realization relevant to materials is that 
of purely attractive pins, similar to those studied in our 
present work. Yet since the properties of elastic lines 
strongly depend on the nature of the pinning centers, 
any conclusions regarding the non-equilibrium relaxation 
properties of disordered type-II superconductors at low 
temperatures that are inferred from models with both 
attractive and repulsive defects should be viewed with 
some scepticism. 



C. Interacting vortex lines without pinning 

In the absence of disorder our system composed of in- 
teracting flux lines evolves toward a regular triangular 
Abrikosov lattice. As we start our simulations by de- 
posing initially straight lines at random positions, large 
displacements of the flux line elements are expected, as 
the system tries to minimize the long-range in-plane re- 
pulsive vortex interaction energy. The ensuing dynamic 
regimes are again nicely captured by the mean-square 
displacement B{t, 0) which takes on values that are two 
orders of magnitude larger than in the absence of inter- 
actions, see the (red) dashed line in Fig. [T] While the 
flux line segments experience these large displacements, 
B(t,0) displays an approximate power law increase with 
time, with an effective exponent of ~ 1.68. Once the 
majority of vortices have reached the vicinity of their fi- 
nal equilibrium positions, the slope of B{t,0) starts to 
gradually decrease. 

The two-time quantities reveal both the initial-time 
regime as well as the crossover at later times, see Fig. |51 
The mean-square displacement s) yields a reason- 
ably good data collapse for the smaller waiting times with 
the scaling exponent —1.68 that follows from the slope 
of B{t,0) in that regime. When the observation time t 
exceeds the crossover time, this scaling breaks down. In- 
stead the growth rate of B{t,s) decreases strongly with 
increasing s, even resulting in a crossing of the curves 
for different waiting times. The behavior of B{t,s) is 
mirrored by that of the vortex density autocorrelation 




FIG. 7: (Color online) Mean-square displacement B{t, 0) vs. 
time t (in MCS) for interacting vortex lines with different val- 
ues of the pinning strength p (system size L = 640; the data 
result from averaging over typically 100 independent runs). 
Due to the long-range interactions, the flux lines aim to max- 
imize their separations, thus yielding values of B{t, 0) that are 
two orders of magnitude larger than in the absence of repulsive 
forces, compare Fig. [J] These displacements are impeded and 
vortex motion eventually stopped by the caging constraints 
of neighboring lines and pinning to attractive point defects. 
The (red) dashed line displays our data for a pure system, 
in the absence of pinning centers. The vertical dashed lines 
indicate the times t = 100, 400, 1600, and 6400 MCS. 



function: For small waiting times s, scaling is achieved 
with exponent 1.68, whereas for larger waiting times the 
decay of the correlation slows down as s increases. From 
these results we infer that the vortex density-density au- 
tocorrelation contains essentially the same physical infor- 
mation as the mean-square displacement. Interestingly, 
the height-height autocorrelation function displays a dif- 
ferent scaling for smaller waiting times, given by the scal- 
ing exponent b = —0.5 of the free line, see Figs. [51; and 
[8ji. This means that during the initial rearrangement of 
the vortex lines the height fluctuations are essentially the 
same as for the free line. Only when the vortices come 
close to their equilibrium positions does the character of 
the correlations change, reflecting the presence of long- 
range repulsive forces. 



D. Interacting vortex lines with pinning 

We are now ready to study the combined effects of 
repulsive flux line interactions and point defect pinning 
during the non-equilibrium relaxation of vortex matter in 
disordered type-II superconductors. We again begin by 
first considering the mean-square displacement B(t^O), 
see Fig. [71 Adding very weak attractive defects, e.g., 
with p — 0.01, has only a very minor effect on the time 
evolution of B{t,0). Strengthening the point pins leads 
to a smaller rate of increase for the mean-square dis- 
placement, see the curve for p = 0.05. At early times the 
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FIG. 8: (Color online) Two-time quantities for interacting 
vortices in the absence of pinning centers: (a) and (b) mean- 
square displacement; (c) and (d) height-height autocorrela- 
tion function; (e) and (f) vortex density-density autocorre- 
lation (with a = 0.05). Data obtained for different waiting 
times s are shown, typically obtained from averaging over 
500 independent simulation runs; the system size is L = 640. 
(a), (c), and (e) display the unsealed data in a log-log plot, 
whereas (b), (d), and (f) show the approximate scaling ob- 
served for not too large waiting times. This scaling regime 
corresponds to the intermediate time window where B{t,0) 
exhibits a power law increase as a function of time, with an 
(effective) exponent ~ 1.68, see Fig. [7] 



flux lines are still displaced from their initial positions, 
as the vortex interactions try to establish an Abrikosov 
lattice. However, at intermediate times these displace- 
ments are impeded by the defects that noticeably slow 
down vortex motion. As a result the system tends to 
a new (quasi-)equilibrium state that balances these two 
competing mechanisms. For even stronger pinning, the 
moving flux lines become rapidly trapped by the disorder 
and the system gradually freezes into a blocked configu- 
ration. In Fig. [7]this is clearly the case for both p = 0.20 
and p = 0.78. 

The most interesting scenario naturally emerges for in- 
termediate defect strengths. Indeed, when the pins are 
very weak, the rearrangement of the flux lines is barely 
affected, and all studied two-time quantities quantita- 
tively display the same behavior as in a pure system. On 
the other hand, when the defects are too strong, the flux 
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FIG. 9: (Color online) Two-time quantities for interacting 
vortices in the presence of attractive pinning centers: (a) 
mean-square displacement; (b) height-height autocorrelation 
function; and (c) vortex density-density autocorrelation (with 
a = 0.05). The pinning strength here is p = 0.05, the system 
size L = 640, and the data result from averaging over typi- 
cally 800 independent simulation runs. The waiting times s 
for the vortex density autocorrelation are the same as those 
for the mean-square displacement. For the height autocor- 
relation function, the waiting times range from s = 25 to 
s = 51200, the values of s being doubled between consecutive 
curves. The (normalized) height-height autocorrelation func- 
tion, after a crossover, displays the same qualitative behavior 
as that encountered in structural glasses. For large values of 
t — s, the height autocorrelation assumes the stretched expo- 
nential form @ with /3 = 0.4, as shown by the dashed line 
overlying the s = 51200 data. 



line elements remain firmly attached to the pinning cen- 
ters, and a frozen configuration ensues. The non-trivial 
behavior at intermediate pinning strengths is studied in 
more detail in Fig. [9] for the case p = 0.05. 

The data for B{t, s) and Cv{t, s) shown in Figs.lH^ and 
IHt are readily understood by comparing them to the cor- 
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responding results for the pure case, see Figs. [8^ and [8^. 
The difference between these data sets is the absence of 
the early-time regime where B{t, 0) has an approximately 
constant slope, see Fig. [T] Consequently the data with 
p = 0.05 do not allow any data collapse, not even for 
the smallest waiting times considered. However, this is 
the only noticeable difference, and the behavior for larger 
waiting times is qualitatively the same as for p = 0, ex- 
cept that the decrease in slope oi B{t,s) for larger values 
of s is stronger when p ^ 0. 

However, a completely different picture emerges for the 
evolution of the normalized height-height autocorrelation 
function. As shown in Fig. [9}3, for waiting times s larger 
than a certain crossover value, C{t, s) exhibits the typical 
two-step relaxation of a structural glass: An initial time- 
translation invariant regime, which corresponds to the 
so-called /3 relaxation in glasses and only depends on the 
elapsed time difference t — s, is followed by a slow decay 
that is usually referred to as a relaxation in the glass 
literature!^ In the long-time limit we can fit this slow 
decay to a stretched exponential 



fir) 



cxp 



(9) 



with T = t ~ s, and a waiting-time dependent decorre- 
lation time td{s). For our different waiting times we ob- 
tain a consistent value f3 « 0.40 for the stretching expo- 
nent in Eq. (|9]). This emergence of a characteristic two- 
step glass-like relaxation is very intriguing. Obviously, 
the flux lines do not settle into a stable microstate even 
after their lateral displacements have become strongly 
reduced owing to the capture by the attractive pinning 
centers and the caging due to their repelling neighboring 
vortices. Instead, as a consequence of the two compet- 
ing relaxation mechanisms, collective dynamics and slow 
decorrelation sets in that yields the typical two-step re- 
laxation dynamics of a glass. 

We also note the intriguing shape of the normalized 
height autocorrelation function in the crossover regime. 
Indeed, at intermediate waiting times C{t,s) displays a 
strongly non-monotonic behavior, with a maximal value 
that even exceeds the value C(s,s) at t = s. This re- 
markable feature points to a fundamental change in the 
nature of the emerging correlations, which is due to the 
trapping of vortex segments in the vicinity of the defects 
and the subsequent balancing of the competition between 
the attractive pinning and the repulsive interactions. 



IV. DISCUSSION AND CONCLUSION 

Our three-dimensional Monte Carlo investigation of re- 
laxation processes in disordered type-II superconductors 
has allowed us to gain a thorough understanding of the 
non-equilibrium properties of these technologically im- 
portant materials. We find the relaxation processes to 
be dominated by the interplay of two competing inter- 
actions, namely the pinning of the flux line elements to 



attractive point defects and the long-range mutual re- 
pulsion of the vortices. This competition generates var- 
ious crossover scenarios that we have discussed system- 
atically. The most interesting regime emerges for pin- 
ning centers of intermediate strength, for which we ob- 
serve a distinguished two-step relaxation and a final slow, 
stretched-exponential decay of the height- height autocor- 
relation function. This behavior is reminiscent of that 
encountered in structural glasses, clearly demonstrating 
that disordered type-II superconductors subject to point 
defects indeed display pronounced glassy behavior at low 
temperatures, again justifying the term "vortex glass" 
for this frustrated pinned low-temperature phase. 

We remark that our results are at variance with re- 
cent investigations based on three-dimensional London- 
Langevin dynamics simulations, where standard aging 
and dynamical scaling behavior of two-time quantities 
was observed, -SLSS, akin to the relaxation features of elas- 
tic lines in a random medium. However, these stud- 
ies, in addition to difference in sample preparation, sys- 
tem size, and boundary conditions in the z-direction, 
employed a coarse-grained continuous random medium 
model of disordered type-II superconductors that be- 
comes adequate near the normal- to superconducting 
transition, but does not realistically capture supercon- 
ducting materials at low temperatures for which isolated 
defects such as oxygen vacancies always induce a local 
suppression of the transition temperature and therefore 
constitute attractive localized pinning centers for vor- 
tices. As our study shows, dynamical scaling no longer 
prevails for purely attractive point pinning centers, but 
instead much richer glassy relaxation dynamics sets in. 

As already mentioned in the Introduction, it is essen- 
tial for investigations of non-equilibrium systems to study 
different dynamics and their algorithmic implementa- 
tions in order to ensure that any ensuing results indeed 
describe actual physical properties of the system rather 
than numerical artifacts. We have therefore recently be- 
gun to implement corresponding London-Langevin dy- 
namics simulations for our elastic line model ([T]) with 
exclusively attractive pinning centers.— Our first tenta- 
tive findings are in complete agreement with the Monte 
Carlo simulation results reported in this paper: They 
too show the emergence of glass-like behavior, with a 
slow, stretched-exponential decay at long times. An in- 
depth analysis of this dynamics is currently in progress; 
this comparative study also aims at matching the differ- 
ent microscopic time scales implicit in Monte Carlo and 
Langevin dynamical simulations. 

Our current study can readily be expanded in various 
directions. Our results are valid in the regime where all 
time-dependent length scales remain small compared to 
the size of the system. However, many transport and 
relaxation experiments are carried out on thin supercon- 
ducting films rather than bulk samples. In our model, 
a finite (small) number of layers introduces a dominant 
new length scale that substantially changes the relaxation 
processes, leading to additional crossover features. We 
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also note that other types of defects can be experimen- 
tally realized, ranging from parallel and splayed colum- 
nar pins to planar defects, and combinations thereof with 
point disorder. It is an open and intriguing problem 
to understand how these different defect configurations 
influence the out-of-equilibrium relaxation processes in 
type-II superconductors. A detailed understanding of the 
relaxation phenomena in superconducting materials may 
facilitate characterization and optimization of samples 
with respect to pinning and flux transport. Finally, in 
all transport applications the flux lines are driven across 
the samples by external currents, which at long times 
yields a non-equilibrium steady state replacing the ther- 
mal equilibrium state that emerges without drive. Fol- 
lowing similar lines as in the present study, one should be 
able to also analyze the relaxation properties of driven 
disordered type-II superconductors in a comprehensive 



manner. We plan to address these and related problems 
in the future. 
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